/**
 * @brief   :
 * @version : 0.1
 * @author  : wanghuan
 * @date    : 2024.11.07
 * @copyright Copyright (c) 2024 by 神工仿真, All Rights Reserved.
 */

/**
 * @brief   : 理论参考对转角场和剪切应变场进行合理插值的厚板元 ，岑（cen）松
 */

#include "TmTPlateCalculator.h"
#include "DataStructure/Input/Property/ShellSecMaterialConstant.h"

using SG::DataStructure::Common::Real;
using SG::Algebra::Matrix;
using SG::FEMSolver::Element::ShellSectionProp1;

void TmTPlateCalculator::SetCoordinate (const Matrix& localCoord,const ShellSectionProp1& secConstant)
{
    auto tp = secConstant.m_t;
    // SecConstant.GetCs(tp, Cs);
    auto factor = secConstant.m_ts * tp;
    auto Cs = Matrix (2, 2, {
        secConstant.m_Cs(0, 0) * factor, secConstant.m_Cs(1, 0) * factor,
        secConstant.m_Cs(1, 0) * factor, secConstant.m_Cs(2, 0) * factor
    });


    // SecConstant.GetCb(tp, Cb);
    const auto& Cb6 = secConstant.m_D;
    Matrix Cb (3,3, {
        Cb6(0,0), Cb6(1,0), Cb6(2,0),
        Cb6(1,0), Cb6(3,0), Cb6(4,0),
        Cb6(2,0), Cb6(4,0), Cb6(5,0)});
    auto IOvr12 = secConstant.FI2 /12.0;
    Cb = Cb * (tp * tp* tp * IOvr12);

    auto b1 = localCoord (1, 1) - localCoord (2, 1);
    auto b2 = localCoord (2, 1) - localCoord (0, 1);
    auto b3 = localCoord (0, 1) - localCoord (1, 1);

    auto c1 = localCoord (2, 0) - localCoord (1, 0);
    auto c2 = localCoord (0, 0) - localCoord (2, 0);
    auto c3 = localCoord (1, 0) - localCoord (0, 0);

    // squre of the length of element side
    std::vector<Real> squreL{ b1 * b1 + c1 * c1, b2 * b2 + c2 * c2, b3 * b3 + c3 * c3 };

    std::vector<Real> delta (3);
    Matrix d(1,2),T1(3,3),T2(2,2);
    for (int iSide = 0; iSide < 3; ++iSide)
    {
        Real lambda;
        if (secConstant.needTrans)
        {
            Real beta;
            if (0 == iSide)
            {
                d(0,0)=localCoord(2,0)-localCoord(1,0);
                d(0,1)=localCoord(2,1)-localCoord(1,1);

                Real dmod=sqrt(d(0,0)*d(0,0)+d(0,1)*d(0,1));
                d(0,0)=d(0,0)/dmod;
                d(0,1)=d(0,1)/dmod;

                beta = acos(d(0,0));
            }
            else if (1 == iSide)
            {
                d(0,0) =-localCoord(2,0);
                d(0,1)= -localCoord(2,1);

                Real dmod=sqrt(d(0,0)*d(0,0)+d(0,1)*d(0,1));
                d(0,0)=d(0,0)/dmod;
                d(0,1)=d(0,1)/dmod;

                beta = acos(d(0,0));
            }
            else if (2 == iSide)
                beta = 0.0;

            Real c = cos(beta); 
            Real s = sin(beta);
            SG::DataStructure::FEM::calcMatCoordTransMatrix(c, s, T1, T2);

            Matrix D_i = T1*Cb*T1.Transpose();
            
            Matrix Cs_i = T2*Cs*T2.Transpose();

            lambda = D_i(1,1)/(Cs_i(1,1)*squreL[iSide]);
        }
        else
            lambda = Cb (0, 0) / Cs (0, 0) / squreL[iSide];

        // end
        delta[iSide] = 6.0 / (1.0 / lambda + 12.0);
    }


    auto A2 = b2 * c3 - b3 * c2;

    // plate data which is independent on the gaussian point
    // 进行了调整（2，3 ）行互换
    m_Bb0 = Matrix (3, 9, { 0, 0,  -b1,  0, 0,  -b2,  0,  0, -b3, 
                            0, c1,   0,  0, c2,   0,  0, c3,  0, 
                            0, b1, -c1,  0, b2, -c2,  0, b3, -c3
                             }) * (-1.0 / A2);

    m_DELTA        = Matrix (3, 3);
    m_DELTA (0, 0) = 1.0 - 2.0 * delta[0];
    m_DELTA (1, 1) = 1.0 - 2.0 * delta[1];
    m_DELTA (2, 2) = 1.0 - 2.0 * delta[2];

    m_DELTA1       = Matrix (3, 3);
    m_DELTA1 (0, 0) = delta[0];
    m_DELTA1 (1, 1) = delta[1];
    m_DELTA1 (2, 2) = delta[2];

    m_oneOverSquareL        = Matrix (3, 3);
    m_oneOverSquareL (0, 0) = 1.0 / squreL[0];
    m_oneOverSquareL (1, 1) = 1.0 / squreL[1];
    m_oneOverSquareL (2, 2) = 1.0 / squreL[2];

    // 变量含义不同， 角度自由度的列互换：如 2和3列互换，并修改第3列的正负号
    m_G = Matrix (3, 9, { 0,   0,   0, -2,  b1, c1,  2, b1, c1, 
                          2,  b2,  c2,  0,   0,  0, -2, b2, c2, 
                         -2,  b3,  c3,  2,  b3, c3,  0,  0,   0,  
                         });
}

Matrix calcFMatrix (Real c1, Real c2, Real c3, Real b1, Real b2, Real b3, const Matrix& oneOverSquareL, Real A2, Real L1, Real L2, Real L3)
{
    auto F = Matrix (3, 3,
                     { c1 * (b2 * L3 + b3 * L2),
                       c2 * (b3 * L1 + b1 * L3),
                       c3 * (b1 * L2 + b2 * L1),
                       -b1 * (c2 * L3 + c3 * L2),
                       -b2 * (c3 * L1 + c1 * L3),
                       -b3 * (c1 * L2 + c2 * L1),
                       (c1 * c2 - b1 * b2) * L3 + (c3 * c1 - b3 * b1) * L2,
                       (c2 * c3 - b2 * b3) * L1 + (c1 * c2 - b1 * b2) * L3,
                       (c3 * c1 - b3 * b1) * L2 + (c2 * c3 - b2 * b3) * L1 });
    return F * oneOverSquareL * (-3 / A2);
}

void TmTPlateCalculator::ComputeB (int iGauss, const Matrix& localCoord, Matrix& Bs, Matrix& Bb, Real& detj) const
{
    Real L2, L3;
    if (iGauss == -1)
    {
        // 计算单元中心（减缩积分）
        L2 = 1.0 / 3;
        L3 = 1.0 / 3;
    }
    else
    {
        const auto& N = m_integForm.GetPoint (static_cast<std::size_t> (iGauss));
        L2            = N.m_x;
        L3            = N.m_y;
    }

    auto L1 = 1 - L2 - L3;

    const auto b1 = localCoord (1, 1) - localCoord (2, 1);
    const auto b2 = localCoord (2, 1) - localCoord (0, 1);
    const auto b3 = localCoord (0, 1) - localCoord (1, 1);
    const auto c1 = localCoord (2, 0) - localCoord (1, 0);
    const auto c2 = localCoord (0, 0) - localCoord (2, 0);
    const auto c3 = localCoord (1, 0) - localCoord (0, 0);

    auto A2 = b2 * c3 - b3 * c2;
    detj    = A2 / 2;

    // bending strain matrix
    auto F = calcFMatrix (c1, c2, c3, b1, b2, b3, m_oneOverSquareL, A2, L1, L2, L3);
    Bb = m_Bb0 + F * m_DELTA * m_G;

    // shear strain matrix
    auto H =
        Matrix (2, 3, { b3 * L2 - b2 * L3, b1 * L3 - b3 * L1, b2 * L1 - b1 * L2,
                        c3 * L2 - c2 * L3, c1 * L3 - c3 * L1, c2 * L1 - c1 * L2 })*(1.0/A2);
    Bs = H * m_DELTA1 * m_G;
}

Matrix TmTPlateCalculator::ComputeStiffness (const Matrix& coord, const Real* const thickness, const SG::FEMSolver::Element::ShellSectionProp1& secConstant) const
{
    Matrix ke (9, 9);
    Matrix Cb;
    Matrix Cs;

    // 弯曲部分
    Matrix Bb;
    Real detj;


    Matrix Bs;
    const auto nGauss = m_integForm.GetPointSize();
    for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
    {
        const auto& tp = thickness[iGauss];
        // SecConstant.GetCb (thickness[iGauss], _OUT Cb);
        const auto& Cb6 = secConstant.m_D;
        Matrix Cb (3,3, {
            Cb6(0,0), Cb6(1,0), Cb6(2,0),
            Cb6(1,0), Cb6(3,0), Cb6(4,0),
            Cb6(2,0), Cb6(4,0), Cb6(5,0)});
        auto IOvr12 = secConstant.FI2 /12.0;
        Cb = Cb * (thickness[0] * thickness[0]* thickness[0] * IOvr12);

        ComputeB (iGauss, coord, _OUT Bs, _OUT Bb, _OUT detj);     
        const Real tmp = detj * m_integForm.GetPoint(iGauss).m_weight; 
        ke += Bb.Transpose() * Cb * Bb * tmp;

        // SecConstant.GetCs (thickness[iGauss], _OUT Cs);
        auto factor = secConstant.m_ts * tp;
        auto Cs = Matrix (2, 2, {
            secConstant.m_Cs(0, 0) * factor, secConstant.m_Cs(1, 0) * factor,
            secConstant.m_Cs(1, 0) * factor, secConstant.m_Cs(2, 0) * factor
        });
        ke += Bs.Transpose() * Cs * Bs * tmp;
    }
    return ke;
}